Introduction

Vegetation physical structure, biodiversity, and topography are interrelated elements controlling key ecosystem properties. Canopy structural complexity—variability in the arrangement of wood and foliar elements in plant canopies (Hardiman et al. 2011, Atkins et al. 2018a)—tends to correlate positively with forest biodiversity (Ehbrecht et al. 2017, Hakkenberg et al. 2018, LaRue et al. 2019) for complex, interrelated reasons. Different tree species may have different growth forms (Verbeeck et al. 2019). More diverse forests, it follows, may have more growth forms represented and thus have more building blocks with which to construct structural complexity (Gough 2020). At the same time, structural complexity also tends to facilitate biodiversity by increasing available niche space for other flora and fauna species (Hansen et al. 1994, Hyde et al. 2006, Torresani et al. 2020).

In the southern Appalachians, complex topography creates steep ecological gradients that drive structural and compositional patterns (Whittaker 1956, Harmon et al. 1984) and produce high levels of biodiversity (Latham and Ricklefs 1993, Gough et al. 2019). Southern Appalachian forests, such as those of Great Smoky Mountains National Park (GRSM), are among the most diverse temperate forests in the world with over 100 woody plant species (Latham and Ricklefs 1993). As such, these forests have long been studied for insights into biodiversity (Whittaker 1956) and related topics including forest productivity (Whittaker 1966, Gough et al. 2019), disturbance (Harmon et al. 1984, Atkins et al. 2020), resource acquisition (Atkins et al. 2018b), habitat fragmentation (Ambrose and Bratton 1990), pollution (Mathias and Thomas 2018), and land cover change (Turner et al. 2003). Perhaps surprisingly, then, we are unaware of earlier studies relating biodiversity and forest structure in this high‐diversity, topographically complex system. Studies conducted elsewhere have found that plant biodiversity tends to decline with increasing elevation (Homeier et al. 2010, Toledo‐Garibaldi and Williams‐Linera 2014), while canopy height and structural complexity also tend to change with elevation (Homeier et al. 2010). For example, along a 3300 m elevational gradient in South America extending from the Amazon lowlands to the tree line in the Peruvian Andes, changes in forest structure were manifested as decreases in canopy height and increases in canopy gap fraction (Asner et al. 2014).

Advances in remote sensing, specifically light detection and ranging (LiDAR), have enabled mapping of forest structure and complexity with unprecedented precision, at plot‐level to regional scales (LaRue et al. 2020). Ecological applications of LiDAR have broadened our understanding of resource acquisition (Stark et al. 2015, Atkins et al. 2018b), allocation strategies (Stovall et al. 2017, 2018a, b, Stovall and Shugart 2018), use efficiencies (Hardiman et al. 2013), drought response (Atkins and Agee 2019, Smith et al. 2019, Stovall et al. 2019, 2020), productivity (Gough et al. 2019), and disturbance history (Fahey et al. 2015, Atkins et al. 2020). However, recent studies of forest structural complexity either focus on broad, landscape‐ to continental‐scale patterns (LaRue et al. 2018, Atkins et al. 2018b, Gough et al. 2019, Fahey et al. 2019), or are limited to characterizations of stand‐scale phenomena (Hardiman et al. 2018, Hickey et al. 2019), without fully considering how forest complexity varies at the stand to regional scale in response to ecological gradients.

Forest structure is most commonly characterized using two‐dimensional indices that describe either vertical or horizontal space only, often ignoring canopy structural traits that describe internal canopy characteristics (Paynter et al. 2018, Fahey et al. 2019) or how forest structure varies simultaneously along both the vertical and horizontal axes. For instance, horizontal structural heterogeneity can be characterized via passive (optical) remote sensing as variation in canopy cover (Morton et al. 2014, Verrelst et al. 2015), while active (radar and LiDAR) remote sensing captures vertical heterogeneity as changes in canopy height (Camaretta et al. 2019). Both metrics of forest structural complexity focus on the outer canopy and ignore internal canopy characteristics (e.g., canopy structural traits such as layering or foliage density). Internal canopy structure can be inferred using statistically derived metrics such as the 90th percentile or standard deviation of LiDAR measured heights above ground level. Statistical metrics derived directly from LiDAR point clouds can be effective and efficient for predictive models of species biodiversity (Goetz et al. 2007), but lack ecological meaning and may be unstable or non‐transferrable across forest types. Recent developments of ecology‐based structural metrics that effectively characterize internal canopy structural traits (Atkins et al. 2018b, Fahey et al. 2019) describe two‐dimensional structure over ground‐based transects and are unable to fully capture three‐dimensional complexity. New metrics need to be developed that synthesize the canopy structural traits such as layering, density, canopy height, and openness in three‐dimensional space to fully characterize structural complexity and its relationship to biodiversity.

To understand relationships between forest structure, biodiversity, and elevation in a topographically complex, high‐biodiversity region, we used three‐dimensional terrestrial LiDAR scanning (TLS) to sample 12 forest plots in Great Smoky Mountains National Park that were previously inventoried by the National Ecological Observatory Network (NEON). We developed two new canopy structural complexity metrics and combined them with traditional TLS‐derived forest structural metrics and vascular plant biodiversity data to address the following questions: (1) How do our newly developed structure metrics compare to traditional metrics? (2) How does forest structure vary with elevation in a high‐biodiversity, high topographic complexity region? (3) How do forest structural metrics vary in the strength of their relationships with vascular plant biodiversity.

Methods

Data collection and post‐processing

Terrestrial LiDAR, or terrestrial laser scanning (TLS), captures detailed 3D information on the forest, enabling the reconstruction and quantitative analysis of forest structure. We visited twelve National Ecological Observatory Network (NEON) distributed vegetation plots (Barnett et al. 2019) located in GRSM (Fig. 1). Plots were chosen based on field accessibility and ancillary NEON data availability at time of scanning. Four TLS scans were acquired in each plot, one at plot center and also at three points around center. Scans were taken between 23 August and 25 August 2018 using two Faro Focus 1203D scanners (FARO, Lake Mary, Florida, USA). Scans were taken at 1/5 resolution and 4 × quality for a total of 28.2 million pulses per scan. All returns with intensity lower than 650 (maximum value = 2100) were removed. Filtering low intensity returns reduces noise and ensures gaps are correctly identified (Stovall et al. 2017).

Filtered scans were used to estimate the vertical distribution of plant material or plant area vegetation density (PAVD) at the plot level. We roughly follow the method described by Calders et al. (2014): (1) elevation normalization, (2) calculate gap probability, and (3) estimate PAVD. The one key difference in the current study is the adoption of a fully 3D digital elevation model to normalize topographic effects, as opposed to simple plane fitting. We found plane fitting was inappropriate in the complex topography of the GRSM national park, resulting in overestimates of total canopy height and a misrepresentation of the vegetation distribution. Topographic normalization followed the following procedure: (1) estimate ground surface, (2) remove spikes and anomalous points, and (3) refit and normalize TLS data with 5 × 5 m resolution surface model. The ground surface is first estimated as the first percentile of height from TLS points in a 5 × 5 m grid over a 200 × 200 m area. The initial model is optimized to remove spikes by excluding all pixels >40% slope. The remaining pixels are refit using a k‐nearest neighbor inverse distance weighting (k = 21), and the resulting surface model is subtracted from the height measurements in the TLS data.

PAVD distributions are derived from a simple calculation of the vertically resolved gap probability (Pgap). Pgap is calculated as

Pgap(θ⎯⎯,𝑧)=1‐∑(𝑧𝑖<𝑧,θ⎯⎯)𝑁(θ⎯⎯)$$\text{Pgap}(\overline{\mathrm{\theta }},z)=1‐\frac{\sum (z_i<z,\overline{\mathrm{\theta }})}{N(\overline{\mathrm{\theta }})}$$

where z is the height above ground and θ is the midpoint of the 5° zenith angle bin used to aggregate the LiDAR returns. The equation essentially calculates the cumulative number of returns per unit normalized height divided by the total number of outgoing laser pulses. In essence, as more vegetation is intercepted, the probability of a laser pulse escaping the canopy decreases.

PAVD is accurately estimated at the hinge angle (a viewing zenith angle of 57.5°) with Pgap. The hinge angle is used to estimate total VAI as this is the angle at which the G‐Function is nearly invariable at 0.5 for all typical leaf angle distributions (MacArthur and Horn 1969, Jupp et al. 2008). A 5° zenith bin between 55° and 60° is used to approximate the hinge angle region. Pgap is converted to the cumulative VAI distribution with respect to height aboveground (L(z)) using:

𝐿(𝑧)≈‐1.1log(Pgap(57.5))$$L(z)\approx ‐1.1\log (\text{Pgap}(57.5))$$

Estimates of L(z) at the top of the canopy approximate the total VAI for a particular sample location. Once cumulative VAI is estimated, Pgap is calculated at 5° zenith angle bins from 0° to 60° and weighted by the sine of the zenith angle to sample the uppermost area of the hemisphere viewable by the laser scanner. Finally, the total VAI in each 1 m vertical bin is estimated as the first derivative of the L(z) curve after weighting with respect to zenith angle.

Metrics

We derived a suite of structural metrics from TLS data to characterize the horizontal and vertical distribution of vegetation (Table 1). As described above, we derived the cumulative VAI and, subsequently, PAVD, for every scan location. Using the PAVD, we calculated foliage density‐dependent median height and 80th and 90th percentiles (i.e., the height at which 80% or 90% of foliage falls below). While these metrics summarize vertical structure, they are unable to quantitatively describe trends in horizontal structure. For this reason, we developed a new metric we call “Depth,” which describes percentiles of canopy penetration at the hinge angle (Fig. 2). Depth percentiles are calculated as the distance from the laser scanner at which a certain percentage of the total number of returns are observed for 10‐degree azimuth bins. A single scan produces 36 estimates of Depth from which a standard deviation can be calculated, which we denote “σDepth” (Fig. 2). For this work, we derived mean and standard deviation values from median Depth estimates for each azimuth angle bin; other percentiles of depth can be calculated (e.g., 80th, 90th), but were colinear with the median estimates, so we use median as the basis for all subsequent analyses. The inherent benefit of using hinge angle‐based Depth percentiles and variability is that this information synthesizes vertical and horizontal structural complexity. Additionally, because Depth and σDepth covary strongly with canopy height (see Results), we also considered both metrics normalized by canopy height, that is, Depth/canopy height and σDepth/canopy height.

Table 1. Definitions of terrestrial LiDAR‐derived structural terms.

Metric	Definition	Source
Vegetation Area Index (VAI)	One‐sided area of all vegetation elements, including leaf, branches, etc	MacArthur and Horn (1969), Jupp et al. (2008), Atkins et al. (2018)
Height Percentiles (e.g., p10, p90)	Percentile height of foliage distribution or Plant Area Vegetation Density (PAVD; e.g., p90 is the height at which 90% of the vegetation is below)	Goetz et al. (2007)
Depth	Percentile of distance or depth specifically at 57.5° zenith angle (e.g., d95 is the distance at which 95% of all LiDAR returns have occurred)	This paper
σDepth	The standard deviation of all mean depth measurements at the 57.5° zenith angle in 10° azimuth bins	This paper

• Depth metrics are introduced in this paper.

Results

Forest structure and complexity

Forest structure varied widely across plots spanning the elevational range of GRSM. Canopy height ranged from ca. 10 to 35 m (Fig. 3). Canopies between 15 and 25 m had the highest variability in PAVD distributions, showing that vegetation layering differs among and within plots (Fig. 3). The high‐elevation locations were dominated by dense conifers with a notable presence of coarse woody debris (e.g., large, downed trees) and standing dead, while the taller locations in lower elevations were more open throughout the canopy (Table 2).

Table 2. Plot elevation, biodiversity, and forest structural metrics.

Plot ID	Elevation (m)	Richness	VAI	p90 (m)	Depth (m)	σDepth (m)
GRSM_001	467	41	3.59	27.25	11.00	6.51
GRSM_002	969	9	3.50	20.67	6.32	2.97
GRSM_006	527	22	3.96	19.25	6.56	3.65
GRSM_013	561	31	3.88	19.75	5.76	3.66
GRSM_014	545	34	3.59	29.75	10.33	5.48
GRSM_016	1778	7	4.02	15.50	4.51	2.76
GRSM_020	608	44	4.02	23.75	7.14	4.19
GRSM_021	587	47	3.25	30.33	11.03	6.37
GRSM_024	624	11	3.54	17.75	3.53	2.99
GRSM_025	1791	9	4.52	11.33	3.34	1.27
GRSM_027	1973	17	3.59	10.00	2.43	1.18
GRSM_029	1786	4	4.02	15.00	5.75	2.04

Notes

• Plot IDs correspond to NEON plot designations. The p90 metric corresponds to canopy height.

Our measures of canopy density (Depth) and structural complexity (σDepth), as well as their counterparts normalized by canopy height, were sensitive to structural variations and effectively summarized variability in canopy structure in the horizontal and vertical dimensions (Fig. 4). Depth summarized information about the density and height of the canopy. Plots with similar Depth could have substantially different σDepth (Fig. 4), as this metric captures how variable Depth is across all 10‐degree azimuth bins. Higher σDepth is indicative of canopies with high variability in vegetation structure: high‐density vegetation and gaps distributed throughout the canopy.

Forest structural metrics considered in this study were moderately to strongly correlated with each other (Fig. 5). Lower stature canopies had higher VAI in GRSM (Fig. 4). VAI was also moderately negatively correlated with both Depth and σDepth, but taller canopies were consistently more structurally complex (Fig. 5). Total VAI was of less relative importance than canopy height in capturing changes in structural variability (Fig. 5). Depth and σDepth were strongly positively correlated with each other (Fig. 5). By design, normalizing Depth and σDepth by canopy height reduced their correlation with canopy height, putatively making them better descriptors of internal canopy structure.

Topographic gradients in the GRSM controlled forest structure and complexity (Fig. 6). Canopy height declined with increasing elevation. VAI was highest in the high‐elevation plots and lowest in the low‐elevation plots. Forest complexity as measured by Depth and σDepth decreased from low to high elevation.

Biodiversity

Plant species richness varied widely among plots, from 4 to 47 species (Table 2). Species richness was not significantly correlated with VAI (Fig. 7; Pearson r = −0.39, P = 0.205), but was strongly positively correlated with canopy height (r = 0.80, P = 0.002), Depth (r = 0.77, P = 0.004), and σDepth (r = 0.84, P < 0.001). Plant species richness was also negatively correlated with elevation (Fig. 8; r = −0.66, P = 0.018). The correlations between biodiversity and forest structural complexity metrics tended to be stronger than the correlation between biodiversity and elevation. Although the strongest correlate of plant species richness was σDepth, the confidence interval on r included the estimates for canopy height, Depth, and σDepth/canopy height.

Discussion

In this study, we developed four new forest structural complexity metrics suited to quantifying aspects of the 3D structure of forests (Depth, σDepth, Depth/canopy height, and σDepth/canopy height) and investigated how they correlate with biodiversity across a high‐biodiversity, topographically complex area in Great Smoky Mountains National Park. Statistically derived LiDAR structural metrics have been limited by either focusing solely on characterizing either horizontal structure (Fischer et al. 2019) or vertical structure (Fotis et al. 2018, Kamoske et al. 2019), or by a lack of true representation of how 3D structure varies with space (Stark et al. 2015, Hardiman et al. 2018, Atkins et al. 2018b). Capturing the complex relationships between biodiversity and physical structure requires the summation of multiple dimensions of variability, including the biological and physical. The bounds of forest structural complexity are likely set by the species pool from which that forest draws. Depth and σDepth metrics simultaneously summarize vegetation structural complexity in both the vertical and horizontal dimensions, describing the internal arrangement of the forest. Our metrics likely correlate strongly with measures of species richness since the greater the number of species, the more likely the species included will be architecturally dissimilar from others.

In this study, we compared our new structural complexity metrics to two of the most widely used metrics of vegetation structure, VAI and canopy height; however, many other vegetation structure metrics have been proposed. Here, we consider how our metrics compare to a selection of these. Ehbrecht et al. (2016, 2017) introduced the effective number of layers (ENL) and stand structure complexity index (SSCI). These metrics are applied to TLS data like our metrics, but represent fundamentally different approaches to quantifying complexity. Our metrics are based on distance from scanner, while the ENL is based on a voxel‐based hit‐grid approach to quantifying vegetation layering. SCCI further incorporates the fractal dimension of polygon slices based on the scanner’s field of view. We expect our metrics to correlate positively with ENL and SCCI and are arguably easier to calculate. Like our metrics, SCCI reflects complexity along both the horizontal and vertical axes. Another metric, foliage height diversity (FHD) was introduced in a classic study of forest complexity (MacArthur and MacArthur 1961) and remains in use today, including as the primary forest complexity product from the GEDI mission (Tang et al. 2019). Similar to ENL, FHD requires discretizing the canopy into height classes and primarily reflects vertical complexity, as opposed to the horizontal complexity emphasized in particular by σDepth.

Potential users should consider the specific aspects of forest structure and complexity represented by our new metrics. Depth simultaneously summarizes information about vegetation density and stature, while σDepth captures internal canopy variability. Since these metrics are more sensitive to internal canopy structure than canopy height, we recommend several TLS acquisitions (>3) per plot to capture variability in complexity—variability itself may also act as a higher order metric of stand heterogeneity. In plots with little spatial variation in stand structure, single scans are likely sufficient to capture the complexity of a plot. Normalizing these new metrics by forest stature enables relative comparisons of stand architecture similar to coefficient of variation (e.g., Are short or tall forest stands relatively more complex?). Further investigations of these metrics are needed in a diverse array of forest types, specifically across gradients of stand age and major environmental factors like temperature, precipitation, soil characteristics, among others.

Across the high‐biodiversity, high‐relief, and topographically complex Great Smoky Mountains, both forest structural complexity and vascular plant biodiversity tended to decline with elevation. In this case, we observed high among‐plot variability in plot‐level species richness (alpha‐diversity), even as high regional species richness (gamma‐diversity) was maintained by substantial species turnover among plots (beta‐diversity). The elevational gradient in GRSM seems to create a microcosm of latitudinal clines in biodiversity, where an increasingly inhospitable environment supports fewer species and constrains plant growth (Whittaker 1956, Hawkins et al. 2003, Gillman et al. 2015). However, inference here is limited, the main caveat being an approximately 1000 m gap in elevation between our low‐ and high‐elevation plots. While many studies of elevation and diversity relationships show monotonic declines in species richness with elevation, there are as many studies showing a humped relationship, with species richness peaking at mid‐elevations (Rahbek 1995, Grytnes and Vetaas 2002). It is important to note that mid‐elevation is relative to the study system—different ranges for what mid‐elevation is would exist when comparing the Appalachians to the Rockies, or either to the Himalayas.

It is unclear whether structural complexity begets biodiversity, or vice versa. Indeed, neither direction of causality strictly excludes the other. Structural complexity could enhance biodiversity by creating niche space that can be filled by other species, or biodiversity could enhance structural complexity due to architectural diversity among species creating more ways to build complexity (Gough et al. 2019, Verbeeck et al. 2019) at higher values of species richness. It may, in principle, be possible to disentangle these effects if a scenario with many plots featuring different subsets of an overlapping species pool. However, with a modest number of plots with high beta‐diversity, this is not possible in our study.

The strong correlations we observed between canopy structural complexity and biodiversity suggest that structural complexity metrics could be used to assay plant biodiversity over large areas. The availability of LiDAR from airborne platforms is increasing due to programs such as NEON (Kao et al. 2012), and new spaceborne sensors such as GEDI (Dubayah et al. 2020) create opportunities for quantifying forest structural complexity at fine spatial resolutions over unprecedented geographic extents. Studies such as ours are important precursors, however, because they explicitly link biodiversity and forest structural complexity, but studies in other systems are needed to understand the generality of forest structure–biodiversity relationship, and determine which metrics, and combinations thereof, are most predictive. The metric foliage height diversity (Tang et al. 2019, Burns et al. 2020), an established GEDI data product, captures vertical structural complexity and is a likely candidate for inclusion in successful predictions of biodiversity from airborne and spaceborne LiDAR sensors.